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Abstract 

We report on a calculation of the quantized energy spectrum and quasiparticle 
eigenfunctions for low lying excitations in the mixed state of clean d-wave 
superconductors. Our study is based on an approximate analytical solution 
of the Bogolubov-de Gennes equations for both rectangular and triangular 
flux lattices with one of the primitive translations chosen parallel to the gap 
node direction. For excitations with momenta close to a certain gap node 
we have obtained a set of eigenfunctions which appear to be extended along 
the chosen gap node direction and localized along the perpendicular one on a 
scale determined by the intervortex distance. The periodic superfluid velocity 
field induces a band structure in the spectrum, which depends essentially on 
the vortex lattice geometry. 

PACS numbers: 74.80.Dm, 74.20.De, 74.60.Ge 
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Recently a great deal of attention has been devoted to the nature of the quasiparticle 
states in isolated vortices and vortex lattices in superconductors with anisotropic pairing. 
This problem is of considerable importance since low-energy quasiparticle excitations in- 
fluence strongly on various static and dynamic properties of the mixed state at low tem- 
peratures. These investigations were stimulated by a large number of experiments which 
provide good evidence for the d-wave symmetry of the order parameter in high-temperature 
superconductors (see [|I| and references therein). For conventional s-wave superconductors 
the low-lying quasiparticle states are bound to the vortex cores as it was first predicted by 
Caroli, de Gennes, and Matricon For an isolated vortex line the eigenvalues for these 
localized states may be written as follows: ~ fiA/{kpC,), where A is the gap value far 
from the vortex axis, ^ is the coherence length at T = 0, kp is the Fermi momentum, and 
the angular momentum quantum number /i is half an odd integer. 

In contrast to conventional superconductors, the density of states (DOS) at low energies 
in d-wave systems is dominated by contributions which come from the regions far from the 
vortex cores and depends essentially on the intervortex distance R^. This unusual fact is 
a direct consequence of the vanishing pair potential in directions of the gap nodes. As a 
result, one obtains the specific magnetic field dependence of DOS at the Fermi level (see 
Refs. |§|-§[), which may be experimentally identified, for instance, in specific heat [^,0 or 
scanning tunneling microscope (STM) |^ measurements. It should be also emphasized that 
there is an important difference between the quantization of the energy spectrum in the 
mixed state of s- and d-wave superconductors: for s-wave systems the energy quantization 
and corresponding localized states exist even in a single isolated vortex line (and they are 
weekly influenced by the presence of neibouring vortices at least for i?„ ^ ^) while for the d- 
wave case the low lying energy spectrum may be quantized only due to the finite intervortex 
distance (see Refs. P|-p!^). The latter conclusion is proved by the numerical solution of the 



Bogolubov-de Gennes (BdG) equations [Tj] which shows that there are no truly localized 
states for a single isolated vortex in a pure d-wave superconductor. 

The most simple theoretical description of the excitation spectrum is based on the semi- 
classical approach which takes account of the Doppler shift of the quasiparticle energy by 
the local superfiuid velocity V^: 



e(k, r) = ±^hW^{k - kpf + A2(k) + h\iY,{r) (1) 

where Vp is the Fermi velocity. Hereafter we assume the Fermi surface to be two- 
dimensional, which is appropriate to high-T^ superconductors, and take the gap function in 
the form = A^kr^ky/kp (the x axis is taken along the [110] crystal direction and thus 
makes an angle 7r/4 with the a axis of the Cu02 planes). We restrict ourselves to the study 
of a vortex lattice under a magnetic field applied parallel to the c axis and, as a consequence, 
the vector V<j is a periodic function of x, y with the periodicity of the vortex lattice. One can 
separate two length scales for quasiparticle wavefunctions: an atomic lengh scale kp^ and a 
characteristic wavelength of a slowly varying envelope /. The semiclassical procedure is in 
general correct when the order parameter and superfiuid velocity are modulated on a scale 
A 3> / . For s-wave superconductors we have / ~ ^ (the minimum spatial extension of wave 
packets made with excitations (|l])) and, as a result, the semiclassical approach fails only for 
the vortex core region. For d-wave systems the I value appears to be angular dependent. In 
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the homogeneous case the spectrum for the low-lying excitations which are close to one of 
the gap nodes (which corresponds, e.g., to the point k = {kp,0)) has the form: 



e{k)c^±Ao^e{k-kfy + k^y/kl, (2) 

where ^ = HVp/Aq. The typical momenta in the x direction are qx ~ 6/{AqC,) and the 
corresponding wavelength is of the order ~ Ao^/e. Comparing the Ix value with the 
intervortex distance (the characteristic length of the superfluid velocity variation), one 
can expect that the semiclassical approach based on the expression (|l|) fails for low energies 
e ^ Aq^/R^. In this case we may assume that the effective potential ^kVs(r) in Eq. (|I]) 
should be averaged over the distances ~ Ix in the x direction (the validity of these quali- 
tative arguments will be proved below). To analyse the problem beyond the semiclassical 
approach one must use the more powerful methods based on either the BdG equations or 
Green's- function techniques (which are equivalent for clean superconductors). The quasi- 
classical limit [kpC, ^ 1) of these theories is known to be represented by the Andreev's 
and Eilenberger equations, respectively. Within these models one should solve the one- 
dimensional quantum mechanical problem for the particle motion along the quasiclassical 
trajectory which is characterized by the impact parameter b = fi/kp and the angle 6 in the 
x — y plane P,P|- P!^JT3 ;T^]. Using the the Bohr-Sommerfeld quantization rule for the angular 



momentum fj,{9, e) one can in principle determine the true quantum levels (see Refs. 0-0) 
For the d-wave case the main difficulty is connected with the description of quasiclassical 
trajectories which make an angle 6 < C,/Rv with the gap node directions. For this angular 
domain the extension of the wave function exceeds the Ry value and quasiparticle states are 
sensitive to the superfluid velocity flelds of all vortices even if the impact parameter is less 
than Ry. The correct description of these trajectories with small 6 values is of consider- 
able importance since it is this angular interval which determines the true quantum levels 
according to the Bohr-Sommerfeld quantization rule. 

It is the purpose of this Letter to report on a calculation of the quantized energy spectrum 
and quasiparticle eigenfunctions for various vortex lattice structures, based on an approx- 
imate analytical solution of the BdG equations for low lying excitations with e ^ AoC,/Rv 
Our approach provides a possibility to take proper account of the influence of the super- 
fluid velocity flelds of all vortices in the vortex lattice on the quasiparticle motion along 
the trajectories almost parallel to gap node directions. We consider the case of pure d- 
wave superconductivity and neglect the effects connected with the possible formation of 
states of mixed symmetry (with coexisting s- and d-wave or (i^2_j^2 and dxy order parameter 
components). The BdG equations can be written as: 

hou{r) + J A(r,r>(r')rfr' = 6u{v) (3) 

- h*ov{r) + J A*(r, r')u{r')dr' = 6v{r) (4) 

where u, v are the particlelike and holelike parts of the quasiparticle wave function. The 
one-particle Hamiltonian ho in the most simple isotropic case takes the form 
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ho = -— V + I— A - E_ 
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where 0o is the flux quantum and Ep is the Fermi energy. The system has been previ- 
ously solved numerically for speciflc lattice models |17,18| and in the continuum limit |l^,|l9|. 
To obtain the analytical solution we follow the well known procedure (see Refs. |T3| , |20[] ) and 
simplify the nonlocal off-diagonal terms using the condition kpC, ^ 1 (or, equivalently, 
Ao ^ Ep). In this case one can search for the solution in the form u = Uexplikpr), 
V = V^ea;p(zki?r), i.e. divide out the fast oscillations on a scale kp^. Then we rewrite A(r, r') 
in terms of the center of mass R = (r + r')/2 and relative coordinates p = r — r' and intro- 
duce the gap function as a Fourier transform with respect to p: A(k, R) = A(^(k)\I^(R). The 
function ^^(R) = fexp{ix) is the d-wave order parameter used in Ginzburg-Landau theory. 
Let us take the two-term expansion for the gap operator: 



J A(r, r')v{r')dr' 



=^^-'vl/A,(k^)r-^^ 



(5) 



where we use the notation {A, B} for the anti-commutator of two operators A and B. In 
order to obtain the Andreev's equations we should keep only the flrst term in this expression. 
Obviously such an approximation is not correct and the second term can not be omitted 
when kp is close to the gap nodes and the flrst term vanishes. Taking, e.g., kpi = {kp,0) 
and introducing a new two-component wave function g = {Uexp{—ix), V) (to eliminate the 
order parameter phase in one obtains the equations linearized in gradient terms: 



\ dx 00 , 




f}9 



hkp Aq 
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where ax,cry,az are the Pauli matrices. We follow here the treatment in Ref. |r3[ and 



neglect the curvature of the Fermi surface. For an isotropic Fermi surface such an approxi- 
mation is valid only for e Aq/ (kpC,) (see Refs. [ pT| , p2| ). However, as mentioned in p!E| , p2| , 
the range of validity of Eq. (|^) may be even larger if the Fermi surface is somewhat flat- 
tened at the nodes. At intermediate magnetic flelds Hd ^ ^ Hc2 (when the intervortex 
distance is much smaller than the London penetration depth) we can assume H = —Hzq 
to be homogeneous and take the gauge A = ifyxo. Here xq, yo, zq are the unit vectors of 
the coordinate system. For a vortex lattice with primitive translations ai,a2 the superfluid 
velocity may be written in the form: 



v.. 



iT^^H ^ [b, Zq] ^^br 



(7) 



where b = nbi + mb2 and bi = 27rif[a2, Zo]/0O) b2 = 27rif[zo, ai]/0o are the primitive 
translations in the reciprocal lattice, n and m run over all possible integers. As we see 
below the solution of Eq. (|^) depends strongly on the flux lattice structure and its orien- 
tation relative to the crystal axes. Contrary to the conventional isotropic superconductors 
(where a hexagonal flux lattice appears to be energetically favourable), for d-wave com- 
pounds previous theoretical works predicted a rich phase diagram, containing triangular, 
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centered rectangular and square lattices with various orientations relative to the ionic lat- 



tice, as a function of magnetic field and temperature |23- Pq] . In principle one should treat 
the problem self-consistently, i.e. calculate at first the energy spectra for various vortex con- 
figurations and then find the lattice structure corresponding to the free energy minimum. 

In this paper we do not solve this self-consistent problem and restrict ourselves to the 
study of energy spectra for several particular lattice structures. Let us choose one of the 
primitive translations (e.g. ai) to be parallel to the gap node direction xq and consider two 
types of lattices: (I) ai = axo, a2 = (^o-yo, Hao? = 0o (rectangular lattice); (II) ai = oxq, 
a2 = a(xo/2 — cryo), Haa^ = 0o (vortices in the unit cell form a shape of an isosceles triangle 
with the base along the x axis). Note that the centered rectangular lattices of the type II 
(with the gradually changing parameter a as a function of H and T values) were found 
to be energetically favourable for a certain region of the H — T phase diagram within a 
generalized London model taking account of nonlocal and nonlinear corrections to the free 
energy PB| , P^ . We also include in our consideration the square lattice tilted by 7r/4 from the 
a axis (type I, a = 1) which is most stable at least for rather high fields and temperatures not 
very close to Tc according to Refs. [p5|-p8| . Such a lattice structure (though elongated in the 



a direction) is close to the one observed experimentally in twinned YBaCuO monocrystals 
by small angle neutron scattering |]2D| and scanning tunneling microscopy 

If we search for the solution of Eq. (|^) as a linear combination of harmonics e'''^^^G{qx, y), 
the periodic functions / and will be responsible for the interaction of harmonics with 
qx and qx + nbix- As we see below, for e < O.SvrAo^/a the wavefunctions G{qx,y) and 
G{qx + nbix, y) do not overlap in the y direction and, consequently, their interaction is 
negligible. In this case the qx momentum component is a good quantum number and one 
can replace the exact periodic potential (|^) in Eq. (D by the effective potential averaged in 
the X direction. The equation for G{qx,y) reads: 

where ^{z) = 2z — (2m + 1) for mKzKm+l, mEZ. We omitted here the 
small corrections of the order ^/a, which are connected with the suppression of the order 
parameter / in vortex cores. It is convinient to introduce F = {ax+d'z)G/2 and dimensionless 



variables y^ H / {acj)^) = z, qx = i^Q^j Ha / (j)Q, e = EAqtt^^J Ha /(f)o. At the m-th interval 
{m < z < m + 1) one obtains: 

dF /I \ 
-tXa,— +(-^z)-E)F + q^axF = 0, (9) 



dz V2 

where = Q — m — 1/2 and A = [nakFC,)'^ is a dimensionless wavelength. The 
procedure analogous to the one used above for the gap node at the point k^^i = {kp, 0) can 
be carried out for the perpendicular gap node direction kp2 = {0,kp). Taking the gauge 



A = —Hxyo and introducing F = (cTx + az)G*/2, —xJ Ha / = z, qy = TiQjH/^a 



e = EA pn^^H/ jac^o), A = ainkp^-' i-2x^Ha/(Po = z, 2g, = nQ^H/{a^o), 2e = 

EAqttC, ^ H/ (o"0o) , A = 4:a{7ikp^)~^ ) for the type I (II) lattices one obtains Eq. (H). The 
symmetry of BdG equations allows us to conclude that the solutions t/3,4, V3,4, £3,4 for two 
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other gap nodes at the points ki^s = {—kp, 0) and kF4, = (0, —kp) can be found by simple 
transformation: 1/3^4 = ¥{2, ^,4 = 25 ^3,4 = ~^i,2- 

It may be useful to note that inside the m-th interval the equations (p|) are equivalent 



to the ones describing the interband tunneling or the one- dimensional (ID) motion of 
a Dirac particle in a uniform electric field and can be solved exactly in terms of parabolic 
cylinder functions. To obtain the energy spectrum one should match these solutions at 
points z = m. This procedure may be essentially simplified using the usual ID quasiclassical 
approach which is valid for A -C 1. Substituting F oc exp{iS{z)) in Eq. (|) we obtain the 
classically allowed (CA) regions (where the function S is real): |$(z)/2 — E\ > \qm\- These 
regions shift in the z direction with a change in the Q value and the quasiparticle motion at 
the m-th interval is classically allowed only for Ig^l < 1/2 + \E\. From the latter condition 
one can see that for rather low energies {\E\ < 1/2) the CA regions corresponding to the 
values Q and Q -|- 2n do not overlap. In this case the interaction of harmonics with and 
+ nbix will be exponentially small due to the exponential decay of the wavefunction in 
classically forbidden regions. Thus, our procedure based on the averaging of the superfiuid 
velocity in the x direction is correct only for low energies \E\ < 1/2, when the solutions are 
essentially localized in the y direction on a scale determined by the intervortex distance. 

Let us now consider the domain \qm\ < 1 and continue with the calculation of quasiclas- 
sical wavefunctions and energy levels corresponding to the quasiparticle motion at the m-th 
interval which contains two CA regions: (i) —l/2<z — m — l/2<E — \qm\ (region A) and 
(ii) E + \qm\ < z — m — 1/2 < 1/2 (region B). Note that the region A (B) exists if the turning 
point zim = m + 1/2 + E — \qm\ (-22m = m + 1/2 + E + \qm\) belongs to the m-th interval. 
The wavefunctions and quantization rules for regions A and B can be written as follows: 

^ ^ Ca,b 




(10) 



1/2±E 



/ ^f^^dt = 7,\{NA,B + lA,B{<im)) (H) 



where t = z — m — 1/2 — E, Na,b is an integer. The ■yA,B values are of the order unity 
and determined by the matching of the expression ([TD|) with the exponentially decaying 
solutions solutions in classically forbidden regions at the (m— l)-th and (m + l)-th intervals. 
This matching procedure results in the following boundary conditions for the wavefunction 
F={F^,F2y. 

, X E-\I2 

= exp{iaA); cosua = -p. — (12) 

Q-m + 1/2 



El 
Fi 

F2 



z=m+l 



E + 1/2 

exp{iaB); cosub = 5777 (13) 

Q — m — 3/2 
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Evaluating the integral in the l.h.s. of Eq. (|T^) one obtains: 




±E 



-ql- qlcosh-' { ' , f'"" ) = 27i\{Na,b + 1a,b) (14) 



The important point is that a set of eigenvalues corresponding to a certain momentum Q 
coincides with that for the momentum Q + 1. Such a periodicity of the energy spectrum is a 
consequence of the periodicity of the potential $(2;) and can be proved exactly from Eq. (^. 
Thus, to analyse the spectrum one can consider only the Q values in the ID Brillouin zone 
— 1/2 < Q < 1/2. In Fig. we display the energy branches (|T^) in the first Brillouin 
zone for the particular case A = 0.01. Each energy branch EA{rn, Na,Q) {EB{m, Nb,Q)) 
considered as a function of the momentum Q has a single minimum (maximum) at the point 
close to the Brillouin zone boundary Q = m + 1/2. Near the first Brillouin zone boundary 
\Q - 1/2| < 1/2 ± E we obtain: 

1 / — 77 7 (g- 1/2)2 J27iX{Na,b + 1a,b) , , 

2 2J27rA(iVA,B + 7AB) 1^ - V^l 



Near the points of intersection of the energy branches (0) in the E — Q plane one should 
take account of the splitting of energy levels resulting from the quasiparticle tunneling 
through classically forbidden regions. As a result we obtain a set of narrow bands separated 
by energy gaps. The effect of tunneling (and, therefore, energy splitting) is exponentially 
small if the characteristic length of the wavefunction decay in a classically forbidden region is 
much smaller than the dimension of this region. In the opposite case the effect of tunneling 
is essential. In particular, we can not neglect the tunneling between the regions A and B for 
Q close to Brillouin zone boundaries, when the distance between the turning points Zim and 
becomes comparable to the characteristic decay length in the region zim < z < Z2m- 
To consider this limit let us take the exact solution at the m-th interval: 



AD_i^_i(r) + 5Z}_i^_i(-r) 
^s%gn{Q -m- l/2)y^(-AD_,^(r) + 5D_,^(-r)) 



(16) 



Here r = ^{z — m — 1/2 — E), /i = g^/(2A), D_j^_i(r) and D^i^^r) are the parabolic 
cylinder functions . Using Eqs. (|T^JT3|jl^) one obtains the system of two equations for A 
and B. The solvability condition for this system results in the quantization rule. Far from 
the turning points the parabolic cylinder functions can be replaced by the corresponding 
asymptotic expressions. Thus, if the turning points zim and Z2m are not too close to the 
m-interval boundaries {zim — m ':$> \qm\, m + 1 — Z2m ^ km I) and the distance Z2m — is 
rather small {\q_m\ ^ V^) then the quantization rule takes the following form: 

Eqn = TTiVA + X{aA - aB)/2 



7 



where N G Z. The spectrum appears to be equidistant only at the Brillouin zone boundaries 
{Q = ±1/2) and the distance between the energy levels depends strongly on the vortex lattice 
geometry. 

Using the solution of Eq. @ given above we obtain four different sets of eigenfunc- 
tions and eigenvalues associated with four gap nodes. The range of validity of our 



approach is restricted by the conditions e < Q.b^Q-n^^ H / (j)Qmin[a and e < 

O.bAQTT^J H/(j)omin[a^^^'^ /2, a^/"^] for type I and II lattices, respectively. In the opposite 



limit the wavefunctions are not localized in the z direction and our approximate solution 
based on the averaging of the superfluid velocity in the gap node direction is no more valid. 
Thus, for quasiparticle states associated with each gap node kp^ there exist two regimes 
with the crossover parameter eRya/i^oO (where R^a is the intervortex distance in the gap 
node direction, a = x,y): (i) for the low energy regime the spectrum has a band structure 
and wavefunctions are extended in the gap node direction and localized in the perpendicu- 
lar one; (ii) for the high energy regime the spectrum is continuous and eigenfunctions are 
extended in both directions. The characteristic energy e of excitations coming into play at 
finite temperatures is of the order ~ T and, consequently, for thermodynamic and trans- 
port quantities one can expect different regimes with the crossover parameters TRyx/i^oO^ 



TRyy/{AQ^). Taking the case R yx ~ Ryy ~ y4^o/H we obtain the crossover parameter 
introduced previously in Refs. 

The quantization of the low energy spectrum should result in the oscillatory behavior 
of DOS as a function of energy with the characteristic energy scales of the order ~ 



/S.QsJhujc/ {aEp) and 8e2 ~ Aq^TiuJcCt / Ep, where ujc = eH/{mc) is the cyclotron frequency. 
Even if we neglect these oscillations and assume to be a continuous variable (more detailed 
study is left as a future problem) the spatial and energy dependence of the local DOS in 
the low energy regime differs qualitatively from the one obtained from the semiclassical 
approach based on the expression (|l|). Taking account of contributions from four gap nodes 
we obtain: 




Nf 

Nrr(e,x,y) = ,, , ^ 

f{z, e) = Mz)/2 + |$(z)/2 + e\ 

where Nj (Nu) is the local DOS for the lattices of the type I (II), Np is 
the density of states at the Fermi level in a normal state. Contour plots of the 
functions ANi{x,y)Np^ J2Hc2/ [t^H) (for the rectangular lattice with a = 2) and 



ANn{x,y)Np^^2Hc2/{-nH) (for the hexagonal lattice with a = ^3/2) are shown in Fig. 
and Fig. |^, respectively. The peaks in the local DOS occur at the points of intersection of 
lines parallel to the x and y axes and passing through the vortex centers. The interesting 
fact is that for type II lattices the peaks in the local DOS appear not only at the vortex 
centers while for rectangular lattices (type I) the coordinates of all peaks coincide with the 
vortex positions. In the vicinity of each vortex center [x, y) the local DOS exhibits a fourfold 
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symmetry which is in a good agreement with the previous studies of a single isolated vor- 
tex in a d-wave superconductor [|l^,|l^ and decreases linearly with an increase of distances 
Iz — and \y — y\ from the vortex center (contrary to the {{x — xY + {y — yY)~^/'^ divergence 
resulting from the semiclassical approach based on Eq. ([|)). The spatially averaged DOS 
has the form 

<'V,>=^,/^|V?+^l|l + -^^| (20) 




If we assume the parameter a to be field independent then the magnetic field dependence 
of the residual DOS (at e = 0) follows the square root behavior which was first predicted in 
Ref . [H] . The deviations from this behaviour may appear if the vortex lattice structure and, 



hence, the a value change with a magnetic field increasing (see Refs. p3| , p^ , p6| -pq] ) and, 
probably, they are the cause for the discrepancies in specific heat measurements 1^,0] • Note 
that these discrepancies may be also explained taking account of the disorder effects @J^. 

In summary, we have described the distinctive features of the low lying quasiparticle 
states in the vortex lattices in d-wave systems. The periodic superfluid velocity field is 
shown to induce the band structure in the low energy part of the spectrum. This band 
structure and corresponding eigenfunctions have been analysed for both rectangular and 
centered rectangular flux lattices tilted by 7r/4 from the a axis. The resulting peculiarities 
of the local DOS are shown to be strongly influenced by the vortex lattice geometry. The 
unusual behaviour of DOS discussed above can be probed by specific heat and STM mea- 
surements and could provide additional arguments in favour of the d-wave pairing in high-T^ 
superconductors. The present study provides a starting point for the analysis of static and 
dynamic properties of the mixed state in various d-wave systems, including, probably, high- 
Tc copper oxides. For this purpose the approach developed above should be generalized 
to describe the quasiparticle states for arbitrary lattice orientations with respect to crystal 
axes. 
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Fundumental Research (Grant No. 97-02-17437). 
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FIGURES 




FIG. 1. The energy branches in the first Brillouin zone for A = 0.01. 
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FIG. 2. Contour plot of the local DOS for the rectangular lattice with a = 2{x' = x{2H / , 
y' = y{2H/^^f-^): (a)^ = 0; (b)^ = 0.27rCAo(2if/</.o)0-^ 
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